                                                                                                                                                                                                                                         % Matlab code to estimate a flexible mixed logit model in wtp space.
% Written by Kenneth Train, first version Oct 21, 2015, 
% Heather Snell added discounting component April 29, 2019.
% Last edits on Jan 4, 2021
 format shortG
 clear all
% ONLY MODIFIED FOR ANA TYPE WITH QH DISCOUNTING AND NO SUBSETS
 
% Do not change the next line. It sets the global variables.
global nset NP XMAT PROBS Z NZ NCS NDRAWS NV BETAS ZTYPE anatype CrossCorr YesGPU COEF NROWS discountfun T modelname subset
% The ANA models are only set up for quasi-hyperbolic type.

% Changing these options may result in the code not working properly.

discountfun = "QuasiHyper";
nset = 2; %two sets of coefficients (one for fully attentive and one for sometimes/always ignored)
subset = 0; %for full data set-not set up for other types
modelname = 'ANA';
anatype = "validation"; %validation has two price coef (used in paper) 
                         %validation1 has one price coef
                         %validation2 has three sets of coef
T=65; %max time periods;

% change diary name and filename as preferred
diary off
delete ANA.out
diary ANA.out
filename = 'ANA.workspace.mat';
CrossCorr=1;
NBins=20;
WantHessian=0;
WantBoot=1;
NReps=250;
NDRAWS=750; %had to slightly reduce with bootstrapping because of memory issues
NGridPts=5001; %had to slighlty reduce with bootstrapping because of memory issues

% ZTYPE=1 for polynomials.
%      =2 for step functions
%      =3 for splines
ZTYPE=3;

PolyOrder=8;
NLevels=6;  
NKnots=8;

R_Range = [0.1765, 0.8825]; %rate
P_Range = [0, 1]; %actually the beta range 
WTP_Range=[-0.2358  0.3460; %att price range
          -0.0311   0.0750; %%ig price range
          -0.3952	1.5; %cap
          -0.5	1.25; %sw
          -0.1659	0.2677;  %att buffer
          -0.25   0.2271;  %ig buffer
          -0.4431   0.4897;  %att infra
          -0.6089   0.2635;  %ig infra
          -0.0501   0.0545; %att jobs
          -0.1647   0.1297; %ig jobs
          -0.2   0.2496; %att quality
          -0.0391   0.15; %ig quality
          -0.75   0.5517; %att wlife
          -1.2372   1.0556];%ig wlife




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(subset==0)
    %load('discounting.data.mat');  %This loads XMAT with the variables below
    % loads different dataset for the ana models
    XMAT = csvread('anadata.csv', 1,0); %skip first row because it has names
    % XMAT is only used for id and csid's
else 
    disp('Not set up for subsets.');
end


FirstRun = 1; %This is used and changed in Bootstapping, but don't change here
Nalts = 3; %number of alts per question
Nquest = 5; %number of questions each person answers
NROWS = size(XMAT, 1);
NP = NROWS/(Nalts*Nquest);
NCS = NP * Nquest; 

% 1. Person number (1-NP)            MUST BE THIS. DO NOT CHANGE.
% 2. Choice situation number (1-NCS) MUST BE THIS. DO NOT CHANGE.
% 3. Chosen alternative (1/0)        MUST BE THIS. DO NOT CHANGE.
% 4. Price 
% 5. CAP
% 6. SW
% 7. Buffer
% 8. Quality
% 9. Jobs
% 10. Infra
% 11. Wlife
% 12. treatment
if(anatype=="elimination")
    NAMES = {'Rate';  'Beta'; 'PriceScale'; 'ASC.CAP'; 'ASC.SW'; 'Buffer';
                'Infra'; 'Jobs'; 'Quality'; 'Wlife'};
elseif(anatype=="validation") %2 price
    NAMES = {'Rate'; 'Beta'; 'AttPrice'; 'IgPrice';'ASC.CAP'; 'ASC.SW'; 'AttBuffer';'IgBuffer';
    'AttInfra'; 'IgInfra'; 'AttJobs'; 'IgJobs'; 'AttQuality'; 'IgQuality'; 'AttWlife'; 'IgWlife'};
elseif(anatype=="validation1") %1 price
    NAMES = {'Rate'; 'Beta'; 'PriceScale'; 'ASC.CAP'; 'ASC.SW'; 'AttBuffer';'IgBuffer';
    'AttInfra'; 'IgInfra'; 'AttJobs'; 'IgJobs'; 'AttQuality'; 'IgQuality'; 'AttWlife'; 'IgWlife'};
elseif(anatype=="validation2") %3 sets
    NAMES = {'Rate'; 'Beta'; 'ASC.CAP'; 'ASC.SW'; 'AttPriceScale'; 'PartPriceSc'; 'IgPriceSc';
             'AttBuffer'; 'PartBuffer'; 'IgBuffer'; 'AttInfra'; 'PartInfra'; 'IgInfra'; 'AttJobs'; 'PartJobs'; 'IgJobs'; 'AttQuality';'PartQuality'; 'IgQuality'; 
             'AttWlife'; 'PartWlife'; 'IgWlife'};
else 
    disp('Not set up for subsets or other ana types.');
end

IDPRICE=4;

NV = size(NAMES, 1);

COEF = [R_Range;P_Range;WTP_Range];
ThisSeed=1234;
rng(ThisSeed);

%Do not change the following lines. They calculate the number of Z variables
if ZTYPE==1; 
    NZ=PolyOrder*NV;
    if CrossCorr==1;
        if discountfun=="QuasiHyper"
            NZ=NZ+(NV-3)*(NV-4)/2; %will not correlate price or rate or beta
        else 
            NZ=NZ+(NV-2)*(NV-3)/2; %will not correlate price or rate
        end
    end
elseif ZTYPE==2
    NZ=(NLevels-1)*NV;
    if CrossCorr==1;
            if discountfun=="QuasiHyper"
               NZ=NZ+2*(NV-3)+(NV-3)*(NV-4)/2; %will not correlate price or rate or beta
            else 
               NZ=NZ+2*(NV-2)+(NV-2)*(NV-3)/2; %will not correlate price or rate
            end
    end
elseif ZTYPE==3
    NZ=(NKnots+1)*NV;
    if CrossCorr==1;
            if discountfun=="QuasiHyper"
               NZ=NZ+2*(NV-3)+(NV-3)*(NV-4)/2; %will not correlate price or rate or beta
            else 
               NZ=NZ+2*(NV-2)+(NV-2)*(NV-3)/2; %will not correlate price or rate
            end
    end
end

%Set the starting values for the coefficients of the Z variables
beta=zeros(NZ,1);  

YesGPU=1;

MAXITERS=5000;

%CONVERENCE TOLERANCE.
% Set =[] to use matlab defaults. default is 1e-6
PARAMTOL=[];
LLTOL=[];
tic;

if(anatype~='validation2')
    CreateDrawsDisWtpProbsQHpref;
    CheckRate;
    if check_ok==1
        CreateZQH;
        disp(' ');
        disp(['Data setup took ' num2str(toc./60) ' minutes.']);
        disp(' ');
        DoEstimationRate;
        paramhat_orig = paramhat;
        if(anatype=="validation")
            Figure_Output16;
        elseif(anatype=="validation1")
            Figure_Output15;
        else 
            Figure_Output10;
        end 
        if WantBoot==1
           Bootana2
        end %end boot
    end %end checkok
else 
    CreateDrawsDisWtpProbsQHpref;
    CheckRate;
    if check_ok==1
        CreateZQH;
        disp(' ');
        disp(['Data setup took ' num2str(toc./60) ' minutes.']);
        disp(' ');
        DoEstimationRate;
        paramhat_orig = paramhat;
        Figure_Output20;
        if WantBoot==1
           Bootana2
        end %end boot
    end %end checkok
    
    
    
end %end if ana type

save(filename)
diary off
